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Abstract 



A Monte Carlo program based on Geant3 has been developed for BESII detector 
simulation. The organization of the program is outlined, and the digitization proce- 
dure for simulating the response of various sub-detectors is described. Comparisons 
with data show that the performance of the program is generally satisfactory. 

Key words: Monte Carlo simulation, Detector response 
PACS: 07.05.Tp 



1 Introduction 



To understand the intrinsic characteristics and performance of a detector, for 
both detector design and physics analysis in experimental high energy particle 
physics, reliable Monte Carlo (MC) simulation is essential. 

The Beijing Spectrometer (BES [1,2]) detector started to take data in 1989 at 
the Beijing Electron Positron Collider (BEPC), which operates in the center- 
of-mass energy range from 2 to 5 GeV. The early simulation program was EGS 
[3] based and thus did not generate hadronic interactions of the secondary 
particles in the detector. Some effects, like fake photons produced by hadrons, 
were not properly described. For the detector response, simple parametric 
models for sub-detectors were used in the simulation. As the program grew, the 
hard-coded structure made it difficult to maintain and make improvements. 

With the upgrade of the BES detector (BESII) [4] and BEPC in 1996, the 
development of a new Monte Carlo program, SIMBES, based on the Geant3 
package [5], was begun. The decision was motivated by many reasons. First, 
larger data samples require better detector simulation to reduce systematic un- 
certainties. Second, some sub-systems of BESII were partly changed and some 
totally new, and it was difficult to make the required software changes within 
the old program structure. Meanwhile, the Geant framework had proven to be 
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a powerful tool for detector description and particle tracking in high energy- 
physics. 

Geant3.21, the last version of Geant3, is used in SIMBES. Much effort is made 
to model the signal generation, since the digitization procedure is highly de- 
tector dependent. BESII data is used for tuning parameters to describe the 
detector response. Many physics channels are compared with data to check the 
simulation, and the agreement is reasonably good. In this paper, the general 
organization and features of SIMBES are described in Section 2. Detailed sim- 
ulations of detector response, as well as comparisons with data, are presented 
in Section 3. 



2 The BESII Simulation Program - SIMBES 

2.1 The BESII Detector 

BESII is a conventional solenoidal magnetic detector that is described in de- 
tail in Ref. [4]. A 12-layer vertex chamber (VC) surrounding the beam pipe 
provides trigger and coordinate information. A 40-layer main drift chamber 
(MDC), located radially outside the VC, provides trajectory and energy loss 
( dE/ dx) information for charged tracks over 85% of the total solid angle. The 
momentum resolution is a p /p = 0.017\/1 + p' 2 (p in GeV/c), and the dE / dx 
resolution for hadron tracks is ~ 8%. An array of 48 scintillation counters 
surrounding the MDC measures the time-of-flight (TOF) of charged tracks 
with a resolution of ~ 200 ps for hadrons. Radially outside the TOF system is 
a 12 radiation length, lead-gas barrel shower counter (SC). This measures the 
energies of electrons and photons over 80% of the total solid angle with an 
energy resolution of a E /E = 22%/\fE (E in GeV). Outside of the solenoidal 
coil, which provides a 0.4 Tesla magnetic field over the tracking volume, is 
an iron flux return that is instrumented with three double layers of counters 
which identify muons with momentum greater than 0.5 GeV/c. 

The geometrical and material descriptions of the detector are given in Section 
3. The experimentally measured field map is used in the simulation. 

2.2 Program Structure 

Built into the design of Geant is the separation of geometry, tracking, and 
digitization. SIMBES respects this feature as illustrated in Fig. 1. In the ini- 
tialization stage, three kinds of data files are loaded: the material and geometry 
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Fig. 1. Schematic structure of SIMBES, dashed boxes mean optional. 



configuration data used to define the detector, the digitization constants to 
describe the detector response, and the realization constants (such as dead 
channels and wire efficiency) to simulate the real detector performance. Input 
values of these parameters can be easily changed by experts or users using 
control cards without relinking the program. 

The BESII event generator package (GENBES) is a stand-alone program that 
can be executed from within SIMBES to provide kinematic quantities for pri- 
mary events. The generators for radiative QED processes [6], selected decays 
of J /if}, ip(2S), and charmed meson pair production [7] were taken from other 
experiments. The BES developed generators, LUND-CRM for inclusive decays 
of J /if) and ip(2S) [8] and PPGEN for some exclusive decays of the if>(2S), as 
well as some standard generators like JETSET, are also included in GENBES. 

At the end of each stage of processing (see Fig. 1), the intermediate output 
(the geometry or hit information in the Geant data structure) can be saved 
to a disk file, and used as input for later SIMBES execution. SIMBES code 
structure is hierarchical and modularized. Therefore, each sub-detector stands 
alone and can be turned on or off by the user. The BESII trigger logic is also 
implemented in a single module; the user has an option to switch it on or off 
when the digitization stage is finished. SIMBES is usually run in batch mode 
to generate raw data and "MC-truth" data. The interactive mode is used for 
graphic event display. 
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2.3 Tracking and Hits 



The tracking of particles through the BESII detector is essentially done ac- 
cording to the Geant standard. Tracking cutoffs and parameters are carefully 
tuned for each sub-detector to provide a compromise between simulation pre- 
cision and computing speed. To correctly simulate electromagnetic shower 
development in the shower counter, the kinetic energy cutoffs of electrons and 
photons have been set to Geant 's limit (10 KeV). 

When a charged particle traverses a sensitive volume (a cell or readout chan- 
nel), the hit information, such as the particle position, energy deposit, etc., 
are saved for the digitization procedure which is discussed in the next section. 



2.4 Hadronic Interaction 



For hadronic interactions, several models have been tried and compared with 
BESII data; each has advantages for specific channels. Fluka [9] and Gcalor 
[10] are the two hadronic models that users can choose in SIMBES. 

Extra photons produced by hadronic interactions in the shower counter are 
studied by a very clean channel ip(2S) — > n + n^J/ip — > 7r + 7r _ / u + /x _ (only 4 
charged tracks in the final state). The neutral tracks observed come mainly 
from pion hadronic interactions in the SC material. Fig. 2 shows comparisons 
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Fig. 2. Number of neutral tracks in the shower counter for 
ip(2S) — ► 7r + 7r - J/V> — > ir + tt~ /i + \i~ decays. Histograms are for MC data 
with (a) Fluka or (b) Gcalor, and points with error bars are for data. 



between data and MC data for different hadronic packages. Gcalor reproduces 
the number of neutral tracks in the SC better than Fluka. 
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In Fluka, inelastic cross sections for anti-nucleon-nucleus annihilation are set to 
zero if the anti-nucleon (n or p) kinetic energy is below 50 MeV (corresponding 
to a momentum of 310 MeV/c); this is, of course, not satisfactory for the 
BEPC energy domain. Therefore, the parametrization given in [11] is used 
instead to approximate the cross-sections for anti-nucleons at low momentum. 
For other particles and for elastic cross-sections, standard Fluka is used. Fig. 
3 shows the energy deposition in the SC for low momentum (less than 0.8 




Fig. 3. Energy deposited by p in the shower counter for J/ip — > ir + ir~pp decays. 
Histograms are for MC data with (a) Fluka or (b) Gcalor, and points with error 
bars are for data. 



GeV/c) anti-protons for J/ip 
than Gcalor's. 



/T + 7r pp decays. Fluka's simulation is better 



3 Simulation of Detector Response 



3. 1 Vertex of event and Vertex Chamber 



The location of the primary event vertex is determined by the interaction 
point (IP) of the beam bunches, which varies with each beam collision. The 
IP location and beam size distributions are obtained from experimental data 
for different running periods, and they are used to sample the primary vertex 
in MC event generation. A T fluctuation caused by the extended beam bunch 
size is taken into account in the TOF simulation (see Section 3.3). 

The vertex chamber (VC), which is built around the beryllium beam-pipe, is 
composed of 12 layers of straw-tubes. Its hit (time) information is used for 
the trigger. When a charged particle traverses a straw-tube, the minimum 
distance between the track and the wire is taken as the drift distance and 
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smeared according to the wire space resolution. The drift time of the ioniza- 
tion electrons is obtained by dividing the drift distance by an effective drift 
velocity. The wire efficiency and the wire sag due to gravity are also taken 
into consideration in this calculation. The VC simulation is checked using the 
decay J/ip — > K S K L [12], which may lose events in the VC trigger because of 
the long decay length of the K$ particles. The agreement between data and 
MC simulation is good when the VC trigger is applied (Fig. 4). 




Fig. 4. Comparison of the Kg (from J/ip — > K$Kl decays) decay length distribu- 
tions between data and MC data. Histograms are for MC, and points with errors 
for data (a) without the trigger and (b) with the trigger. 



3.2 Main Drift Chamber 



The main drift chamber, which is used to determine trajectories and measure 
energy losses of charged particles, is the main tracking detector of BES11. It 
is cylindrical with 10 super layers with 4 sense wires (axial or stereo) in each 
superlayer. The axial layer is described by Geant's TUBE volume. For stereo 
layers, the shape of the hyperbolic tube (HYPE) in Geant cannot reproduce 
the geometry. Instead, the subtraction technique is used on two HYPEs for the 
stereo-layer description. In the simulation, the wires (sense, field and guard) 
are not put into SIMBES individually, but an equivalent material inside the 
MDC, taken as a mixture of gas and wires, is used. 

The energy loss calculated by Geant for a thin layer is not used in simulating 
the energy loss by charged tracks in the MDC. Instead, the energy loss (Q) at 
a sense wire is simulated according to the distribution obtained from real data, 
smeared with an experimentally determined Landau distribution. Fig. 5 shows 
dE / dx distributions for different particles, where dE / dx is the truncated 
mean of the energy loss for each track. 
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3Cj[ %k 

Fig. 5. Top: dE/dx (pulse-height) versus momentum for (a) MC and (b) data. 
Bottom: The Xi (( dE/ dx) measured - ( dEj dx) expected ) / a) distributions for (c) ir 
(d) K particles from inclusive J/ijj decays. Histograms are for MC data, and points 
with error bars are for data. 

Using wire hit information, the minimum distance between a track and each 
wire is calculated including the effects of the Lorentz angle (a^ = 26°) and 
gravitational sagging of the wires. The drift distance of the ionization is taken 
as the minimum distance and is smeared by the wire space resolution (a w ). It 
is observed from data that the residual of the drift distance is dependent on 
Q, and even for a given Q, it cannot be simply fitted with a single Gaussian 
(see Fig. 6). Therefore, the drift distance is smeared according to a combined 
distribution of two Gaussians (fit from data) with the amount of smearing 
determined by the Q value and tracking layer. The wire efficiency is included 
in the a w simulation. When there are multiple hits in a cell, the time signal is 
taken as the minimum drift time, while the pulse height is determined from 
the sum of energy losses. 

Fig. 7 shows momentum distributions (after track reconstruction) for data and 
MC data. The comparison of momentum resolutions for data and Monte Carlo 
is shown in Table 1. The pion tracking efficiency (see Fig. 8) for simulated 
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Fig. 6. Distributions of drift distance residuals for different pulse heights, Q, 
corresponding to ADC counts of (a) 100-500, (b) 900-1100, (c) 1900-2100, and (d) 
3700-4000. Points with error bars are from Bhabha events at the J/ip, and curves 
are fits with combined distributions of two Gaussians. 




P GeV/c P GeV/c 

Fig. 7. Momentum distributions for (a) protons and (b) anti-protons from J/ift — > pp 
decay. Histograms are for MC, and points with error bars are for data. 

data is also found to be consistent with data. In general, an error of 2% per 
track is taken as the systematic error from the tracking efficiency. Alternative 
algorithms for the MDC wire resolution simulation have also been tested, and 
their performances are slightly worse; however, they can be used to determine 
the systematic uncertainties from the MDC wire resolution simulation. 
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Table 1 



Channel 


omc (MeV/c) 


a DT (MeV/c) 


e + e~ — ► at J/ij; 


47.01 ± 0.06 


46.70 ± 0.08 


J/ip ^ pp 


33.22 ± 0.28 


34.18 ± 0.12 
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Fig. 8. 7r tracking efficiency in the MDC versus momentum. The first three points 
are from J ftp — > AA — > 7r + 7r _ pp, and the others are from J/-0 — > 7r + 7r~7r . 



3.3 Time of Flight System 



The TOF system provides precise time measurements for particle identifica- 
tion and the trigger. It consists of two parts, the barrel (BTOF) and the 
end-caps (ETOF). The BTOF, located just outside of the MDC, employs 48 
scintillation counters with photo-multiplier tubes (PMT) at each end. Each 
ETOF contains 24 scintillators with PMTs at one end. In SIMBES, the TOF 
system is described simply by trapezoidal scintillator bars with small gaps 
between adjacent ones. 

After the track passes through the scintillator, the TOF hit information con- 
tains the particle's position, Z, at the entry point, its energy deposition, Ei oss , 
in the scintillator, and the time of flight, Tf, from the IP. Ei oss is determined 
using energy loss with 5-ray generation and applying Birk's law for saturation, 
while for Tf, the time spread due to the sizable beam bunch is taken into ac- 
count. A parameterized model is used to simulate the ADC and TDC outputs; 
it assumes the light produced at the entry point propagates with attenuation 
to the two ends. The amplitude and arrival time are obtained according to 

Q = / x E loss x exp(Z/L atte ),T = 7) + Z/v c + T C (Q, Z), (1) 



where f is a conversion constant, L atte the attenuation length (about 3m), v c 
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the effective light velocity in scintillator, the last term in T is a correction 
for Q and Z which is obtained from the calibration. The time walk effect is 
included in T c . 

ADC and TDC counts are generated by smearing Q and T according to their 
resolutions, and a Q threshold is applied to the signal. The resolutions for 
each counter are obtained using Bhabha events; in the smearing, the time 
resolution dependence with amplitude and hit position is taken into account. 
For multiple hits in a single counter, taking the time walk effect as a good 
approximation, only the hits with Q over threshold contribute to the TDC 
count, while all pulse heights are summed to give the ADC count. 

Fig. 9 shows the AT distributions for Bhabha and dimuon events, and the 




Fig. 9. AT (T measure d— T expec ted) distributions for (a) Bhabha and (b) dimuon events 
at the J/tp energy. Histograms are for MC data, and points with error bars are for 
data. 

pion particle identification (PID) efficiency for the BTOF as a function of 
momentum is shown in Fig. 10. The simulation is in a good agreement with 
data, and an uncertainty of 1% per track is usually taken as the systematic 
error for PID with the TOF. 



3.4 Shower Counter 



The shower counter consists of two parts, the barrel (BSC) and the end- 
caps (ESC). In this section, BSC is taken as an example to demonstrate the 
simulation procedure. The cylindrical BSC consists of 24 layers of gas tubes 
interleaved with 23 layers of lead absorbers. There are 560 cells (tubes made 
of Aluminum) in each gas layer. In order to reduce the number of electronic 
readout channels, layers with the same angle are grouped into six readout 
layers in the r direction. To save computation time, the absorber layer is 
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Fig. 10. Pion identification efficiency as a function of momentum for the TOF from 
J/ij) vr+vr-vr ). 

modeled using a mixture of Al-Pb-Al sandwiches as its material, while the 
sensitive layer (cell) is filled with gas working in self-quenching streamer (SQS) 
mode. The support ribs and some insensitive regions around them are properly 
described in SIMBES. 

When a charged particle hits a gas tube, the energy deposited (total charge) 
is sampled according to the SQS spectrum which can be parameterized by a 
Landau distribution. Corrections are made to this charge for the angle depen- 
dence of the track and for different layers. The pulse heights at the two ends 
of the tube are derived by charge division, and the final ADC outputs are 
obtained by smearing the pulse heights with a Gaussian which accounts for 
the contribution from the electronics. In the case of multiple hits in a tube, if 
two hits are close enough, they are merged as one hit. 

The simulated hit profiles (number of hits in different layers) of an electro- 
magnetic shower are in close agreement with data. Fig. 11 shows the energy 
deposition (using Gcalor) in the BSC for different particles after the energy 
calibration with Bhabha events. Fig. 12 [13] shows the ratio of photon detec- 
tion efficiency in the BSC for data and MC data. The comparison indicates 
reasonable agreement, and 2% per photon is taken as the systematic error for 
the photon detection efficiency. 

3.5 Muon Counter 

The muon identifier is the outermost component of BES and is composed 
of three layers of iron absorber and three layers of proportional chambers. 
There are 189 chambers distributed in an octagonal structure, each chamber 
consists of 8 proportional gas tubes which are arranged into two sub-layers. 
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and dimuon events at the J/ip, (c) and (d) for 7r and K from inclusive J/V> decays, 
and (e) and (f) for p and p from J/ip — > pp. Histograms are for MC data, and points 
with error bars are for data. 

The complex geometrical structure and magnetic field [14] of the muon system 
are well described in SIMBES. 



In the simulation, the hit position (Z) along the tube is smeared by the res- 
olution (obtained from data), the ADC signals at the two ends are obtained 
from charge division of the energy deposit from Geant, and the wire efficiency 
from data is also included. Fig. 13 shows the comparison of the hit positions 
in the muon chamber for data and MC data. The agreement is satisfactory. 



4 Summary 



A Monte Carlo program SIMBES has been developed for BESII simulation. 
Special efforts are made to accurately simulate the response of the detector. 
For each sub-detector, some key physics quantities are compared with BE- 
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Fig. 12. The ratio of photon detection efficiency between data and MC (cdt/cmc) 
versus photon energy, eor and emc are calculated using J/ip — ► 7r + 7r _ 7r° decays. 
The solid line is a fit to the points. 
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Fig. 13. Hit position residuals for three Muon layers; the top, middle, and bottom 

plots correspond to layers 1, 2, and 3, respectively. Muons are from e + e~ — ► 

at J/ip events. Histograms are for MC data, and points with error bars are for data. 

SII data to test the simulation quality. It is shown from many checks that the 
overall performance of the current SIMBES is satisfactory and that the results 
from simulation are generally consistent with data. On the other hand, there 
is still room for further improvement. For example, the digitization of MDC 
should be regarded as an 'effective' simulation where many physics mecha- 
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nisms are hidden, and noise in the MDC is not considered at present. 

The new Monte Carlo package (SIMBES) has significant impact on physics 
results. As an example, the branching ratio of J/ip — > 7r + 7r~7r° measured by 
BES is (2.10 ± 0.12)% [15], which is significantly higher than the PDG value 
of (1.50 ± 0.20)% [16]. In this measurement, it is found that the efficiency 
using a two-Gaussian resolution function in the MDC wire simulation is about 
20% lower than that using a single Gaussian function, and that the efficiency 
including hadronic interactions is about 10% lower than without them. 
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